home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YDBLVDP.C < prev    next >
C/C++ Source or Header  |  1988-11-28  |  2KB  |  88 lines
  1. /********************** YDblVdp for double vanderPol *************************/
  2. /********************* (C) 1986,7,8 by JAMES A. YORKE ************************/
  3.  
  4. #include "yinclud.h"
  5.  
  6. DEDblVdp(k, Y)            /*  2 coupled vanderpol oscillators -- see also
  7.                    map 23 */
  8. double  k[],
  9.         Y[];
  10. {
  11.     int     i;
  12.     double  tval,        /* t is a global variable */
  13.             s,
  14.             X1,
  15.             X2,
  16.             Y1,
  17.             Y2;
  18.     int     base;
  19.     tval = Y[0];
  20.     X1 = Y[1];
  21.     Y1 = Y[2];
  22.     X2 = Y[3];
  23.     Y2 = Y[4];
  24.     s = rho * sin(C4 * tval);
  25.     k[0] = 1;
  26.     k[1] = Y1;
  27.     k[2] = C1 * Y1 * (1 - X1 * X1) - C2 * X1 - C3 * X1 * X1 * X1 + s + rho * X2;
  28.  
  29.     k[3] = Y2;
  30.     k[4] = C5 * Y2 * (1 - X2 * X2) - C6 * X2 - C7 * X2 * X2 * X2 + s + rho * X1;
  31.  
  32.     if(num_lyap > 0) {
  33.         for(i = 0; i < num_lyap; i++) {
  34.             base = lyapzero + vec_dim * i;
  35.             k[base] = Y[base + 1];
  36.             k[base + 1] = (-2 * C1 * X1 * Y1 - C2 - 3 * C3 * X1 * X1) * Y[base]
  37.                 + C1 * (1.- X1 * X1) * Y[base + 1]
  38.                 + rho * Y[base + 2];
  39.             k[base + 2] = Y[base + 3];
  40.             k[base + 3] =
  41.                 rho * Y[base]
  42.                 + (-2 * C5 * X2 * Y2 - C6 - 3 * C7 * X2 * X2) * Y[base + 2]
  43.                 + C5 * (1 - X2 * X2) * Y[base + 3];
  44.         }
  45.     }
  46. }
  47. v2_mod() {
  48.     double  moduloAB();
  49.     y[0] = moduloAB(y[0], 0., twopi / C4);/* time */
  50.     modFlag = NO;
  51.  /* bring y[1] and y[3] into -pi,pi etc. */
  52.     y[1] = moduloAB(y[1], -pi, pi);
  53.     y[3] = moduloAB(y[3], -pi, pi);
  54. }
  55.  
  56.  
  57. initDblVanderpol() {
  58.     extern int      DEDblVdp(), v2_mod ();
  59.     zeroth = 1;
  60.     vec_dim = 4;        /* the dimension of the Lyapunov vectors =
  61.                    phase space dim */
  62.     num_lyap = 0;        /* the number of Lyapunov numbers to be
  63.                    computed <= vec_dim */
  64.     lyapzero = 5;        /* y[lyapzero] is the zeroth coord of the
  65.                    zeroth lyapunov vector */
  66.     dim = lyapzero + num_lyap * vec_dim;
  67.  /* needed for rungekutta */
  68.     X_upper = 2.2;        /* x scale */
  69.     X_lower = -2.2;
  70.     Y_upper = 3.3;        /* y scale */
  71.     Y_lower = -3.3;
  72.     X_coord = 1;
  73.     Y_coord = 2;
  74.     C1 = 2./ 5.;
  75.     C2 = 0.;
  76.     C3 = 1.;
  77.     C4 =.6;
  78.     C5 = 1./ 5.;
  79.     C6 = 0.;
  80.     C7 = 1.;
  81.     rho = 1.0;
  82.     steps_per_cycle = 40;
  83.     step = (twopi / C4) / steps_per_cycle;
  84.     DE = DEDblVdp;        /*  DE is a pointer to a function   */
  85.     modPointer = v2_mod;
  86. }
  87.  
  88.